In [1]:
using Plots
using Random, Distributions
using LinearAlgebra
using StaticArrays
using SpecialFunctions
using BenchmarkTools
using ProgressMeter
const Vec3 = SVector{3, Float64}
Out[1]:
Lo primero, es pensar cómo representar las partículas y los nodos del arbol del algoritmo de Barnes Hut. Lo ideal sería crear un tipo de datos que represente a cada uno y guarde la información relevante. Julia no soporta las clases, por lo tanto, vamos a crear estructuras de datos para representar partículas y nodos.
In [2]:
mutable struct Body
id::Int
force::SVector{3, Float64}
position::SVector{3, Float64}
velocity::SVector{3, Float64}
mass::Float64
Body() = new()
Body(id, force, position, velocity, mass) = new(id, force, position, velocity, mass)
end
Body(id, position, velocity, mass) = Body(id,SVector(0.0,0.0,0.0),position,velocity,mass)
function reset_force(body::Body)
body.force = SVector(0.0, 0.0, 0.0)
end
mutable struct Cell
cell_opened::Bool # If true it means it has children
children::Array{Cell, 1}
body_present::Bool
body::Body
mass::Float64
cm::SVector{3, Float64}
min_bounds::SVector{3, Float64}
max_bounds::SVector{3, Float64}
end
Cell()=Cell(false,Array{Cell,1}(undef,8),false,Body(),0.0,SVector(0.0,0.0,0.0),SVector(0.0,0.0,0.0),SVector(0.0,0.0,0.0))
Out[2]:
Ahora, vamos a crear el esquema de integración temporal para el movimiento de cada partícula. Se utiliza PEFRL, que es un integrador simplectico de orden 4 en la posición. También existe una variación de este algoritmo que se llama VEFRL que es de orden 4 en la velocidad.
In [3]:
function Bmove_r(body::Body,dt::Float64,cte::Float64=1)
body.position+=body.velocity*dt*cte
end
function Bmove_v(body::Body,dt::Float64,cte::Float64)
body.velocity+=body.force*dt*cte
end
function step(particles::Array{Body,1},Δt::Float64,θ::Float64) #PEFRL
const1 = 0.1786178958448091 #ζ
const3 = -0.6626458266981849e-1 #χ
const4 = -0.2123418310626054 #λ
const2 = (1-2*const4)/2 #(1-2λ)/2
const5 = 1-2*(const3+const1) #1-2*(ζ+χ)
#Update tree
(universe_min, universe_max)=get_universe_size(particles)
root = Cell()
root.min_bounds = universe_min
root.max_bounds = universe_max
for particle in particles
insert_body(root,particle)
Bmove_r(particle,Δt,const1)
end
#Time step (Los for son mas rapidos que el broadcast)
for particle in particles
reset_force(particle)
compute_force(root, particle, θ)
end
for particle in particles
Bmove_v(particle,Δt,const2)
Bmove_r(particle,Δt,const3)
end
for particle in particles
reset_force(particle)
compute_force(root, particle, θ)
end
for particle in particles
Bmove_v(particle,Δt,const4)
Bmove_r(particle,Δt,const5)
end
for particle in particles
reset_force(particle)
compute_force(root, particle, θ)
end
for particle in particles
Bmove_v(particle,Δt,const4)
Bmove_r(particle,Δt,const3)
end
for particle in particles
reset_force(particle)
compute_force(root, particle, θ)
end
for particle in particles
Bmove_v(particle,Δt,const2)
Bmove_r(particle,Δt,const1)
end
end
function simulate(particles::Array{Body,1}, steps::Int64; Δt::Float64 = 0.1, θ::Float64 = 0.5, Nframes=200, show_progress=false)
# Array to store simulation data
frames = Array{Body,1}[]
draw=round(Int32,steps/Nframes)
p = Progress(steps, enabled=show_progress, showspeed=true)
# Store initial conditions
push!(frames, deepcopy(particles))
for i = 1:steps
#Perform time step integration
step(particles, Δt, θ)
#Save simulation data every draw steps
if i%draw==0
push!(frames, deepcopy(particles))
end
next!(p)
end
return frames
end
Out[3]:
Lo siguiente es hacer las funciones que definen el árbol y la lógica necesaria para crearlo y añadir las partículas
In [4]:
function cell_contains_position(cell::Cell, position::SVector{3, Float64})
if position[1] < cell.min_bounds[1] || position[1] > cell.max_bounds[1]
return false
end
if position[2] < cell.min_bounds[2] || position[2] > cell.max_bounds[2]
return false
end
if position[3] < cell.min_bounds[3] || position[3] > cell.max_bounds[3]
return false
end
true
end
function insert_body(cell::Cell, body::Body)
if !cell.body_present && !cell.cell_opened
cell.body = body
cell.mass = body.mass
cell.cm = copy(body.position)
cell.body_present = true
else
if !cell.cell_opened
half_sides = (cell.max_bounds .- cell.min_bounds) ./ 2
for i in 1:8
child = Cell()
cell.children[i] = child
shifts = Vec3([((i - 1) >> (c - 1)) & 1 for c in 1:3])
child.min_bounds = cell.min_bounds .+ (shifts .* half_sides)
child.max_bounds = cell.max_bounds .- ((1 .- shifts) .* half_sides)
# Insert old body
if (cell.body_present && cell_contains_position(child, cell.body.position))
insert_body(child, cell.body)
cell.body_present = false
end
end
cell.cell_opened = true
end
# Insert new body
for child in cell.children
if cell_contains_position(child, body.position)
insert_body(child, body)
break
end
end
# Update mass and center of mass of the cell
cell.cm = (cell.mass .* cell.cm .+ body.mass .* body.position) ./ (cell.mass + body.mass)
cell.mass += body.mass
end
end
function get_universe_size(bodies::Array{Body})::Tuple{Vec3, Vec3}
universe_min = [Inf, Inf, Inf]
universe_max = [-Inf, -Inf, -Inf]
for body in bodies
for c in 1:3
if body.position[c] < universe_min[c]
universe_min[c] = body.position[c]
end
if body.position[c] > universe_max[c]
universe_max[c] = body.position[c]
end
end
end
(Vec3(universe_min), Vec3(universe_max))
end
Out[4]:
Ahora, hacemos las funciones que calculan la fuerza sobre cada molécula teniendo en cuenta el árbol.
In [5]:
function add_force(body::Body, position::SVector{3, Float64}, mass::Float64; ϵ::Float64 = 0.002)::Nothing
r = position-body.position
dist = norm(r)+ϵ
dist_cubed = dist * dist * dist
f = G*mass/dist_cubed #f/m aceleration
body.force += f*r
return nothing
end
function compute_force(cell::Cell, body::Body, θ=0.5)::Nothing
cell_body::Body = cell.body
if (cell.body_present && cell_body.id === body.id) || cell.mass === 0.0
return
end
if !cell.cell_opened
add_force(body, cell_body.position, cell_body.mass)
else
size = sum(cell.max_bounds .- cell.min_bounds)
d = norm(cell.cm - body.position)
if size / d < θ
# Treat cell as a body
add_force(body, cell.cm, cell.mass)
else
for child in cell.children
compute_force(child::Cell, body::Body)
end
end
end
end
Out[5]:
Ahora, quisieramos poder graficar las partículas, los octantes y además hacer una animación de la simulación, por lo tanto definimos las siguientes funciones
In [6]:
function show_particles(particles, θ::Float64 = 30.0, ϕ::Float64 = 30.0; Lim=1)
scatter([Tuple(p.position) for p = particles],
lims = (-abs(Lim), Lim),
camera = (θ, ϕ),
size = (500, 500),
label = "",
background_color = :black,
marker = (:circle, 3, 0.8, :white, stroke(0)))
end
function animate_frames(frames; frame_rate::Int64 = 1, file_name = "animation", θ::Float64 = 30.0, ϕ::Float64 = 30.0, quadview::Bool = false, Lim=1)
anim = Animation()
for i = 1:frame_rate:length(frames)
if quadview
show_quadview(frames[i], Lim=Lim)
else
show_particles(frames[i], θ, ϕ, Lim=Lim)
end
frame(anim)
end
gif(anim, "$file_name.gif")
end
function show_quadview(particles; Lim=1)
scatter([Tuple(p.position) for p = particles],
lims = (-abs(Lim), Lim),
camera = (30, 30),
size = (800, 800),
label = "",
background_color = :black,
marker = (:circle, 3, 0.8, :white, stroke(0)),
layout = 4,
subplot = 1)
scatter!([Tuple(p.position) for p = particles],
lims = (-abs(Lim), Lim),
camera = (0, 90),
label = "",
marker = (:circle, 3, 0.8, :white, stroke(0)),
subplot = 2)
scatter!([Tuple(p.position) for p = particles],
lims = (-abs(Lim), Lim),
camera = (90, 0),
label = "",
marker = (:circle, 3, 0.8, :white, stroke(0)),
subplot = 3)
scatter!([Tuple(p.position) for p = particles],
lims = (-abs(Lim), Lim),
camera = (60, 30),
label = "",
marker = (:circle, 3, 0.8, :white, stroke(0)),
subplot = 4)
end
Out[6]:
Ya teneiendo todo lo necesario para la simulación, viene la parte física. La función que calcula la aceleración y la que inicializa el sistema. Es en estas dos funciones es donde se guarda la información de las unidades y la física del problema.
1. Modifique el código de tal manera que se utilice un sistema de unidades adecuado para la descripción de un sistema galáctico. Sugiero un sistema con [M⊙ ], [kpc] y [Gyr].¶
In [7]:
using Unitful, UnitfulAstro
using PhysicalConstants
G = PhysicalConstants.CODATA2018.G
Out[7]:
In [8]:
G = uconvert.(u"kpc^3 *Msun^-1 *Gyr^-2",G)
Out[8]:
In [9]:
G = ustrip(G)
Out[9]:
2. Modifique el codigo para permitir que cada particula pueda tener una masa diferente (puede ser asignada aleatoriamente, dentro de un rango de masas determinado, en el proceso de inicialización ).¶
El struct con el que se inicializan las partículas permite colocar una masa diferente a cada una, por ejemplo, veamos como se ve una inicialización con un número aleatorio de partículas
In [10]:
Random.seed!(123)
function rand_particles(num_particles::Int64, Black_Hole_Mass, R) #2D galaxy
particles = Body[]
u1 = Uniform(0,1)
u50 = Uniform(1,50)
u11 = Uniform(-1,1)
push!(particles, Body(1,[0.0, 0.0, 0.0], [0.0, 0.0, 0.0], Black_Hole_Mass)) #super massive black hole
for i = 2:num_particles
θ = 2π*rand(u1) #Uniform random numbers between [0,1]
r = R*rand(u1) #Uniform random numbers between [0,R]
z = R*rand(u11)/4
rr = SVector(r*cos(θ), r*sin(θ), z)
v = cross(SVector(0,0,sqrt(G*Black_Hole_Mass/(norm(rr)^3))), rr) #v=Ω×r
mass = rand(u50) #Uniform random numbers between [1,50]
push!(particles, Body(i, rr, v, mass))
end
particles
end
bodies = rand_particles(800, 4e6, 30)
display(show_particles(bodies,Lim=30))
3. Agregue al código una función que permita obtener una gráfica de velocidad orbital de todas las partículas vs. distancia radial al centro de la distribución.¶
In [35]:
function graph_vel(particles; Lim=40)
data_r = zeros(length(particles))
data_v = zeros(length(particles))
for i in 1:length(particles)
data_r[i]= norm(particles[i].position)
data_v[i]= norm(particles[i].velocity)
end
plot = scatter(data_r,data_v)
xlabel!("r")
ylabel!("v")
xlims!((0,Lim))
return plot
end
Out[35]:
4. Modelo 1. Agregue al código una opción que permita un modelo que incluya un agujero negro supermasivo con masa M en el centro de la galaxia y una distribución aleatoria de partículas con velocidades angulares iniciales Keplerianas,¶
\begin{equation} \Omega_k(r)=\sqrt{\frac{GM}{r^3}} \end{equation}5. Modelo 2. Agregue al código una opción que permita un modelo en el que las partículas se distribuyan con un perfil de densidad superficial exponencial¶
\begin{equation} \Sigma(r)=\Sigma_0e^{-r/R_d} \end{equation}donde la longitud de escala puede estimarse del orden de $R_d = 0,1R$ con R el radio de la galaxia y donde las velocidades tangenciales inciales estén dadas por \begin{equation} v_0^2=4\pi G\Sigma_0R_dy^2(I_0(y)K_0(y)-I_1(y)K_1(y)) \end{equation}
con $y = \frac{r}{2R_d}$ y donde $I_n$ y $K_n$ son las funciones modificadas de Bessel (estas pueden ser implementadas utilizando librerías).
Después de realizar estas modificaciones, utilice los modeles 1 y 2 implementados para ilustrar el comportamiento de la Vía Láctea utilizando alrededor de 5000 estrellas con masas entre 1M⊙ y 50M⊙ , en un intervalo de tiempo de 13 Gyr y utilizando los siguientes parámetros
Radio de la Galaxia: 30 kpc
Masa del agujero negro central: 4 × 10^6 M⊙
En el caso del Modelo 2 puede escalar la distribución de masa como considere adecuado, siempre y cuando tenga razones físicas para sustentar las suposiciones que realice. Obtenga también la curva de velocidad de las estrellas en el instante final de la evolución para cada uno de los modelos.
Para poder distribuir las masas cumpliendo con el perfil exponencial, sabemos que si x distribuye uniforme entre a y b, entonces el valor esperado es $\mu=\frac{a+b}{2}$. Entonces, la masa total de la galaxia va a ser en promedio $M_T=N\frac{a+b}{2}$ donde N es el número de estrellas y a y b el rango de posibles valores de las masas de las estrellas. Además, la integral sobre la distribución de densidad en la galaxia debe ser igual a la masa total, entonces \begin{split} \begin{equation} M_T=\int_0^{2\pi}\int_0^R\Sigma_0e^{-r/R_d}drd\theta\\ M_T=2\pi\Sigma_0(R_d-R_de^{-R/R_d})\\ \Sigma_0=\frac{M_T}{2\pi(R_d-R_de^{-R/R_d})} \end{equation} \end{split}
Para poder distribuir las partículas, estas se van colocando por anillos para que la densidad tenga la masa adecuada, de modo que la masa en un anillo de la galaxia que va de $r_1$ hasta $r_2$ es \begin{split} \begin{equation} M=\int_0^{2\pi}\int_{r_1}^{r_2}\Sigma_0e^{-r/R_d}drd\theta\\ M=2\pi\Sigma_0R_d(e^{-r_1/R_d}-e^{-r_2/R_d})\\ \end{equation} \end{split} Entonces, por cada región de la galaxia, se bebe asegurar que haya la cantidad correcta de masa.
In [12]:
t = 13 #Gyr
N = 2000 #particles
R = 30 #Galaxy radii kpc
M = 4e6 #Black hole mass
Out[12]:
In [13]:
function rand_kepler_particles(num_particles::Int64, Black_Hole_Mass, R) #2D galaxy
particles = Body[]
u1 = Uniform(0,1)
u50 = Uniform(1,50)
push!(particles, Body(1,SVector(0.0, 0.0, 0.0), SVector(0.0, 0.0, 0.0), Black_Hole_Mass)) #super massive black hole
for i = 2:num_particles
θ = 2π*rand(u1) #Uniform random numbers between [0,1]
r = R*rand(u1) #Uniform random numbers between [0,R]
z = 0
rr = SVector(r*cos(θ), r*sin(θ), z)
v = cross(SVector(0,0,sqrt(G*Black_Hole_Mass/(norm(rr)^3))), rr) #v=Ω×r
mass = rand(u50) #Uniform random numbers between [1,50]
push!(particles, Body(i, rr, v, mass))
end
particles
end
function rand_exp_particles(num_particles::Int64,R; partitions=10) #2D galaxy
particles = Body[]
Rd = 0.1*R
u1 = Uniform(0,1)
u50 = Uniform(1,50)
Mt = num_particles*(50-1)/2
Σ0=Mt/(2*pi*(Rd-exp(-R/Rd)))
dr = R/partitions
r = 0
for i in 1:partitions
u = Uniform(r,r+dr)
M = 2*pi*Σ0*Rd*(exp(-r/Rd)-exp(-(r+dr)/Rd))
while(M>0)
θ = 2π*rand(u1)
m = rand(u50)
rr= rand(u)
y = rr/(2*Rd)
v0 = sqrt(4*pi*G*Σ0*Rd*y^2*(besseli(0,y)*besselk(0,y)-besseli(1,y)*besselk(1,y)))
push!(particles, Body(i,SVector(r*cos(θ), r*sin(θ), 0.0), SVector(v0*sin(θ), -v0*cos(θ), 0.0), m))
M-=m
end
r+=dr
end
println("Total Number of particles: ", length(particles))
particles
end
Out[13]:
Definamos los parámetros para la simulación
In [14]:
Lim1 = 35
particles1 = rand_kepler_particles(N,M,R)
show_particles(particles1,Lim=Lim1)
Out[14]:
Ahora, hagamos las simulaciones
In [15]:
dt = 0.005 #Gaños
n = round(Int64, t/dt)
@time frames1 = simulate(particles1, n, Δt=dt, Nframes=300, show_progress=true)
println("Model 1 done.")
In [16]:
animate_frames(frames1, frame_rate=5, quadview = true, file_name = "animation1", Lim=Lim1)
Out[16]:
In [37]:
graph_vel(frames1[length(frames1)],Lim=35)
Out[37]:
In [18]:
Lim2 = 35
particles2 = rand_exp_particles(N,R,partitions=N/2)
show_particles(particles2,0.0,90.0,Lim=Lim2)
Out[18]:
In [30]:
t = 90 #Gyr
dt = 0.01 #Gyr
n = round(Int64, t/dt)
@time frames2 = simulate(particles2, n, Δt=dt, Nframes=300, show_progress=true)
println("Model 2 done.")
In [31]:
animate_frames(frames2, frame_rate=5, quadview = true, file_name = "animation2", Lim=Lim2)
Out[31]:
In [32]:
graph_vel(frames2[length(frames2)],Lim=45)
Out[32]:
Comparación con el algoritmo tradicional¶
Para ver si enserio se gana tiempo con este algoritmo hagamos la comparación con el algoritmo tradicional. Para ello, vamos a usar la implementación del ejercicio 3.
In [22]:
mutable struct Molecule
#Position
position::SVector{3, Float64}
#Velocity
velocity::SVector{3, Float64}
#Mass
mass::Float64
#force
force::SVector{3, Float64}
end
function clear_f(body::Molecule)
body.force=SVector(0.0, 0.0, 0.0)
end
function move_r(body::Molecule,dt::Float64,cte::Float64=1)
body.position+=body.velocity*dt*cte
end
function move_v(body::Molecule,dt::Float64,cte::Float64=1)
body.velocity+=body.force*dt*cte/body.mass
end
function calc_f(body::Array{Molecule,1}; ϵ::Float64=0.02)
for i in 1:length(body)
for j in i+1:length(body)
dr = body[j].position-body[i].position
d = norm(dr)
F = -G*body[i].mass*body[j].mass*dr/((d+ϵ)^3)
body[i].force+=-F #3 Newton Law
body[j].force+=F
end
end
end
function step_PEFRL(body::Array{Molecule,1},Δt::Float64)
const1 = 0.1786178958448091 #ζ
const3 = -0.6626458266981849e-1 #χ
const4 = -0.2123418310626054 #λ
const2 = (1-2*const4)/2 #(1-2λ)/2
const5 = 1-2*(const3+const1) #1-2*(ζ+χ)
for particle in body
move_r(particle,Δt,const1)
clear_f(particle)
end
calc_f(body)
for particle in body
move_v(particle,Δt,const2)
move_r(particle,Δt,const3)
clear_f(particle)
end
calc_f(body)
for particle in body
move_v(particle,Δt,const4)
move_r(particle,Δt,const5)
clear_f(particle)
end
calc_f(body)
for particle in body
move_v(particle,Δt,const4)
move_r(particle,Δt,const3)
clear_f(particle)
end
calc_f(body)
for particle in body
move_v(particle,Δt,const2)
move_r(particle,Δt,const1)
end
end
function Propagate(particles::Array{Molecule,1}, steps::Int64; Δt::Float64 = 0.1, Nframes=200, show_progress=false)
frames = Array{Molecule,1}[]
draw=round(Int32,steps/Nframes)
push!(frames, deepcopy(particles)) #save initial condition
p = Progress(steps, enabled=show_progress, showspeed=true)
for i in 1:steps
#Perform time step integration
step_PEFRL(particles,dt)
#Save simulation data every draw steps
if i%draw==0
push!(frames, deepcopy(particles))
end
next!(p)
end
return frames
end
function rand_kepler_particles2(num_particles::Int64, Black_Hole_Mass, R) #2D galaxy
u1 = Uniform(0,1)
u50 = Uniform(1,50)
particles = Molecule[]
push!(particles, Molecule([0.0, 0.0, 0.0], [0.0, 0.0, 0.0], Black_Hole_Mass, [0.0,0.0,0.0])) #super massive black hole
for i = 2:num_particles
θ = 2π*rand(u1) #Uniform random numbers between [0,1]
r = R*rand(u1) #Uniform random numbers between [0,R]
z = 0
rr = [r*cos(θ), r*sin(θ), z]
v = cross([0,0,sqrt(G*Black_Hole_Mass/(norm(rr)^3))], rr) #v=Ω×r
mass = rand(u50) #Uniform random numbers between [1,50]
push!(particles, Molecule(rr, v, mass, [0.0,0.0,0.0]))
end
particles
end
Out[22]:
In [23]:
particles3 = rand_kepler_particles2(50,4e6,30)
show_particles(particles3,Lim=31)
Out[23]:
In [24]:
@time frames3 = Propagate(particles3, 2500, Δt=0.005, Nframes=300, show_progress=true)
println("Brutte force done")
In [25]:
animate_frames(frames3, frame_rate=5, quadview = true, file_name = "animation3", Lim=31)
Out[25]:
Hagamos una comparación con los tiempos de ejecución.
In [26]:
N = [3, 9, 20, 60, 100, 150, 200, 400, 700, 1000, 2000, 5000]
timeB1 = []
timeB2 = []
timeB3 = []
timeN = []
for i in N
particlesB = rand_kepler_particles(i,4e6,30)
particlesN = rand_kepler_particles2(i,4e6,30)
tB1 = @elapsed frames = simulate(deepcopy(particlesB), 10, Δt=0.01, θ=0.1, Nframes=2)
tB2 = @elapsed frames = simulate(deepcopy(particlesB), 10, Δt=0.01, θ=0.4, Nframes=2)
tB3 = @elapsed frames = simulate(deepcopy(particlesB), 10, Δt=0.01, θ=0.6, Nframes=2)
tN = @elapsed frames = Propagate(deepcopy(particlesN), 10, Δt=0.01, Nframes=2)
push!(timeB1, tB1)
push!(timeB2, tB2)
push!(timeB3, tB3)
push!(timeN, tN)
end
In [27]:
plot(N,timeB1,label="Barnes-Hut; θ=0.1")
plot!(N,timeB2,label="Barnes-Hut; θ=0.4")
plot!(N,timeB3,label="Barnes-Hut; θ=0.6")
plot!(N,timeN,label="Brute force")
xlabel!("N")
ylabel!("Time [s]")
Out[27]:
In [28]:
plot(N,timeB1,label="Barnes-Hut; θ=0.1")
plot!(N,timeB2,label="Barnes-Hut; θ=0.4")
plot!(N,timeB3,label="Barnes-Hut; θ=0.6")
xlabel!("N")
ylabel!("Time [s]")
Out[28]:
In [29]:
plot(N,timeB1,label="Barnes-Hut; θ=0.1",yaxis=:log,xaxis=:log)
plot!(N,timeB2,label="Barnes-Hut; θ=0.4",yaxis=:log,xaxis=:log)
plot!(N,timeB3,label="Barnes-Hut; θ=0.6",yaxis=:log,xaxis=:log)
plot!(N,timeN,label="Brute force",yaxis=:log,xaxis=:log)
xlabel!("N")
ylabel!("Time [s]")
Out[29]:
In [ ]: